Multistate Bennett Acceptance Ratio Research Articles

Bayesian Multistate Bennett Acceptance Ratio Methods.

The multistate Bennett acceptance ratio (MBAR) method is a prevalent approach for computing the free energies of thermodynamic states. In this work, we introduce BayesMBAR, a Bayesian generalization of the MBAR method. By integration of configurations sampled from thermodynamic states with a prior distribution, BayesMBAR computes a posterior distribution of free energies. Using the posterior distribution, we derive free energy estimations and compute their associated uncertainties. Notably, when a uniform prior distribution is used, BayesMBAR recovers the MBAR's result but provides more accurate uncertainty estimates. Additionally, when prior knowledge about free energies is available, BayesMBAR can incorporate this information into the estimation procedure by using nonuniform prior distributions. As an example, we show that by incorporating the prior knowledge about the smoothness of free energy surfaces, BayesMBAR provides more accurate estimates than the MBAR method. Given MBAR's widespread use in free energy calculations, we anticipate BayesMBAR to be an essential tool in various applications of free energy calculations.

Improved Sampling of Adaptive Path Collective Variables by Stabilized Extended-System Dynamics.

Because of the complicated multistep nature of many biocatalytic reactions, an a priori definition of reaction coordinates is difficult. Therefore, we apply enhanced sampling algorithms along with adaptive path collective variables (PCVs), which converge to the minimum free energy path (MFEP) during the simulation. We show how PCVs can be combined with the highly efficient well-tempered metadynamics extended-system adaptive biasing force (WTM-eABF) hybrid sampling algorithm, offering dramatically increased sampling efficiency due to its fast adaptation to path updates. For this purpose, we address discontinuities of PCVs that can arise due to path shortcutting or path updates with a novel stabilization algorithm for extended-system methods. In addition, we show how the convergence of simulations can be further accelerated by utilizing the multistate Bennett's acceptance ratio (MBAR) estimator. These methods are applied to the first step of the enzymatic reaction mechanism of pseudouridine synthases, where the ability of path WTM-eABF to efficiently explore intricate molecular transitions is demonstrated.

Revisiting the pseudo-supercritical path method: An improved formulation for the alchemical calculation of solid-liquid coexistence.

Alchemical free energy calculations via molecular dynamics have been applied to obtain thermodynamic properties related to solid-liquid equilibrium conditions, such as melting points. In recent years, the pseudo-supercritical path (PSCP) method has proved to be an important approach to melting point prediction due to its flexibility and applicability. In the present work, we propose improvements to the PSCP alchemical cycle to make it more compact and efficient through a concerted evaluation of different potential energies. The multistate Bennett acceptance ratio (MBAR) estimator was applied at all stages of the new cycle to provide greater accuracy and uniformity, which is essential concerning uncertainty calculations. In particular, for the multistate expansion stage from solid to liquid, we employed the MBAR estimator with a reduced energy function that allows affine transformations of coordinates. Free energy and mean derivative profiles were calculated at different cycle stages for argon, triazole, propenal, and the ionic liquid 1-ethyl-3-methyl-imidazolium hexafluorophosphate. Comparisons showed a better performance of the proposed method than the original PSCP cycle for systems with higher complexity, especially the ionic liquid. A detailed study of the expansion stage revealed that remapping the centers of mass of the molecules or ions is preferable to remapping the coordinates of each atom, yielding better overlap between adjacent states and improving the accuracy of the methodology.

Understanding the sources of error in MBAR through asymptotic analysis.

Many sampling strategies commonly used in molecular dynamics, such as umbrella sampling and alchemical free energy methods, involve sampling from multiple states. The Multistate Bennett Acceptance Ratio (MBAR) formalism is a widely used way of recombining the resulting data. However, the error of the MBAR estimator is not well-understood: previous error analyses of MBAR assumed independent samples. In this work, we derive a central limit theorem for MBAR estimates in the presence of correlated data, further justifying the use of MBAR in practical applications. Moreover, our central limit theorem yields an estimate of the error that can be decomposed into contributions from the individual Markov chains used to sample the states. This gives additional insight into how sampling in each state affects the overall error. We demonstrate our error estimator on an umbrella sampling calculation of the free energy of isomerization of the alanine dipeptide and an alchemical calculation of the hydration free energy of methane. Our numerical results demonstrate that the time required for the Markov chain to decorrelate in individual states can contribute considerably to the total MBAR error, highlighting the importance of accurately addressing the effect of sample correlation.

Alchemical Free Energy and Hamiltonian Replica Exchange Molecular Dynamics to Compute Hydrofluorocarbon Isotherms in Imidazolium-Based Ionic Liquids.

Ionic liquids (ILs) have shown promise for applications that leverage differential gas solubility in an IL solvent, e.g., gas separations. Although most available literature provides Henry's law constants, the ability to efficiently estimate full isotherms is important for engineering design calculations. Molecular simulation can be used as a tool to predict full isotherms of gas in ILs. However, particle insertions or deletions in a charge-dense IL medium and the sluggish conformational dynamics of ILs present two sampling challenges for these systems. We therefore devised a method that uses Hamiltonian replica exchange (HREX) molecular dynamics (MD) combined with alchemical free energy calculations to compute full solubility isotherms of two different hydrofluorocarbons (HFCs) in imidazolium-based IL binary mixtures. This workflow is significantly faster than the Gibbs ensemble Monte Carlo (GEMC) simulations which fail to deal with the slow conformational relaxation caused by the sluggish dynamics of ILs. Multiple free energy estimators, including thermodynamic integration, free energy perturbation, and multistate Bennett acceptance ratio method, provided consistent results. Overall, the simulated Henry's law constant, isotherm curvature, and solubility trends match experimental results reasonably well. We close by calculating the full solubility isotherms of two HFCs in IL mixtures that have not been reported in the literature, demonstrating the potential of this method to be used for solubility prediction and setting the stage for future computational screening studies that search for the "best" IL to separate azeotropic HFC mixtures.

The reconciliation between the experimental and calculated octanol-water partition coefficient of 1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine using atomistic molecular dynamics: an open question

The octanol-water partition coefficient of 1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine (DPPC) was investigated using atomistic molecular dynamics simulations via thermodynamic integration and multistate Bennett acceptance ratio methods. The GAFF and CHARMM36 force fields were used with six water models widely used in molecular dynamics simulations. The OPC4 water model provided the best agreement with the experimental octanol-water partition coefficient of DPPC using the two force fields. However, there is still plenty of room for improvement in water models with correct estimation of surface tension that uses better and suitable non-bonded interaction parameters between water-water and water-DPPC. The Gibbs free energy of transferring DPPC from octanol to water phase was calculated to be 19.8 ± 0.3 and 20.2 ± 0.3 kcal mol−1, giving a partition coefficient of 14.5 ± 0.4 and 14.8 ± 0.3 for the GAFF and CHARMM36 force fields, respectively. This study reinforces the importance of developing new water models that reproduce experimental surface tensions to reconcile the water-water and water-DPPC non-bonded interactions and the existing discrepancy between experimental measurements of amphiphilic molecules that are important in many areas of scientific applications and industry such as biophysics, surfactant, colloids, membranes, medicine, nanotechnology, and food and pharmaceutical industries, and so on. It raises two important open questions: Is the experimental octanol-water partition coefficient of DPPC reliable? Or is its calculation accurate using the OPC4 water model? With respect to the experimental measurements, there may be non-treated aspects such as the formation of aggregates in aqueous phase and limit of detection of the applied method. And, in the calculation, some effects are not possible to be considered in a correct way or viable time such as calculating quantum effects, sampling all conformations, considering phase transitions, and correctly evaluating the intermolecular forces to estimate an accurate surface tension. Communicated by Ramaswamy H. Sarma

Stochastic Approximation to MBAR and TRAM: Batchwise Free Energy Estimation

The dynamics of molecules are governed by rare event transitions between long-lived (metastable) states. To explore these transitions efficiently, many enhanced sampling protocols have been introduced that involve using simulations with biases or changed temperatures. Two established statistically optimal estimators for obtaining unbiased equilibrium properties from such simulations are the multistate Bennett acceptance ratio (MBAR) and the transition-based reweighting analysis method (TRAM). Both MBAR and TRAM are solved iteratively and can suffer from long convergence times. Here, we introduce stochastic approximators (SA) for both estimators, resulting in SAMBAR and SATRAM, which are shown to converge faster than their deterministic counterparts, without significant accuracy loss. Both methods are demonstrated on different molecular systems.

Iterative Landmark-Based Umbrella Sampling (ILBUS) Protocol for Sampling of Conformational Space of Biomolecules.

Computer simulations of biomolecules such as molecular dynamics often suffer from insufficient sampling. Due to limited computational resources, insufficient sampling prevents obtaining proper equilibrium distributions of observed properties. To deal with this problem, we proposed a simulation protocol for efficient resampling of collected off-equilibrium trajectories. These trajectories are utilized for the initial mapping of the conformational space, which is later properly resampled by the introduced Iterative Landmark-Based Umbrella Sampling (ILBUS) method. Reconstruction of static equilibrium properties is achieved by the multistate Bennett acceptance ratio (MBAR) method, which enables efficient use of simulated data. The ILBUS protocol is geometry-based and does not demand any additional collective variable or a dimensional-reduction technique. The only requirement is a set of suitably spaced reference conformations, which serve as landmarks in the mapped conformational space. Additionally, the ILBUS protocol encompasses an iterative process that optimizes the force constant used in the umbrella sampling simulation. Such tuning is an inherent feature of the protocol and does not need to be performed by the user in advance. Furthermore, even the simulations with suboptimal force constants can be used in estimates by MBAR. We demonstrate the feasibility and the performance of this approach in the study of the conformational landscape of the alanine dipeptide, met-enkephalin, and adenylate kinase.

Absolute binding free energies of mucroporin and its analog mucroporin-M1 with the heptad repeat 1 domain and RNA-dependent RNA polymerase of SARS-CoV-2

The peptide Mucroporin and its analog Mucroporin-M1 were studied using the molecular docking and molecular dynamics simulation of their complexation with two protein targets, the Heptad Repeat 1 (HR1) domain and RNA-dependent RNA polymerase (RdRp) of SARS-CoV-2. The molecular docking of the peptide-protein complexes was performed using the glowworm swarm optimization algorithm. The lowest energy poses were submitted to molecular dynamics simulation. Then, the binding free energies of Mucroporin and its analog Mucroporin-M1 with these two protein targets were calculated using the Multistate Bennett Acceptance Ratio (MBAR) method. It was verified that the peptides/HR1 domain complex showed stability in the interaction site determined by molecular docking. It was also found that Mucroporin-M1 has a much higher affinity than Mucroporin to the HR1 protein target. The peptides showed similar stability and affinity at the NTP binding site in the RdRp protein. Additional experimental studies are needed to confirm the antiviral activity of Mucroporin-M1 and a possible mechanism of action against SARS-CoV-2. However, here we indicate that Mucroporin-M1 may have potential antiviral activity against the HR1 domain with the possibility for further peptide optimization. Communicated by Ramaswamy H. Sarma

Computational analysis of the tryptophan cation radical energetics in peroxidase Compound I

Three well-characterized heme peroxidases (cytochrome c peroxidase = CCP, ascorbate peroxidase = APX, and Leishmania major peroxidase = LMP) all have a Trp residue tucked under the heme stacked against the proximal His heme ligand. The reaction of peroxidases with H2O2 to give Compound I results in the oxidation of this Trp to a cationic radical in CCP and LMP but not in APX. Considerable experimental data indicate that the local electrostatic environment controls whether this Trp or the porphyrin is oxidized in Compound I. Attempts have been made to place the differences between these peroxidases on a quantitative basis using computational methods. These efforts have been somewhat limited by the approximations required owing to the computational cost of using fully solvated atomistic models with well-developed forcefields. This now has changed with available GPU computing power and the associated development of software. Here we employ thermodynamic integration and multistate Bennett acceptance ratio methods to help fine-tune our understanding on the energetic differences in Trp radical stabilization in all three peroxidases. These results indicate that the local solvent structure near the redox active Trp plays a significant role in stabilization of the cationic Trp radical.Graphical abstract

Thermodynamic and structural description of relative solubility of the flavonoid rutin by DFT calculations and molecular dynamics simulations

Understanding the solubility of flavonoids is of relevance in view of the broad application such as antioxidant and antitumor agents. The solubility of flavonoid rutin has been experimentally investigated in a homologous series of alcohols, acetone, n-heptane and water solvents. In this work theoretical evaluation of Gibbs free energy of solvation was performed at the Density Functional Theory (DFT) and Molecular Dynamics (MD) levels, with the ωB97x-D functional and multistate Bennett acceptance ratio (MBAR) method, respectively. The degree of solubility of a given solute among various solvents is dictated by solute-solvent interactions, which may be represented by the solvation energy. Our calculated relative energy values due to solvent effects, using quantum chemical methods (PCM and SMD solvent models) do not follow strictly the observed solubility trend. However, the PCM model including solute-solvent dispersion interaction, solute-solvent repulsion interaction and solute cavitation energy contribution to the total energy (named DIS,REP,CAV) yielded an overall agreement with experimental solubility profile, with exception of ethanol solvent. A similar satisfactory accordance with experiment was found for MBresults, with a discrepancy also found for 1-butanol, which deserve further investigation. Analysis of spatial and radial distribution functions (SDF and RDF) for pure solvents and the analysis of the contributions of vdW and electrostatic terms for the total free energy reveal that a balance between dispersive and electrostatic interactions between the rutin and the solvent is necessary to solubilize the flavonoid, which corroborates with the values of these energy contributions to the solvation-free energy calculated by the MBAR method. The use of explicit solvent molecules in DFT calculations is precluded due to the size, flexibility, and many sites for solvent interactions of rutin molecule. Our results indicate that in general DFT-PCM (DIS,REP,CAV) and MBAR levels of calculations can satisfactorily predict the solvation free energy and solubility trend for flavonoid rutin along a series of solvents, and also provide a reasonable estimate of relative solubilities of other flavonoids, information that can be relevant for further molecular modeling studies.

Computational studies of protein–drug binding affinity changes upon mutations in the drug target

AbstractMutations that lead to drug resistance limit the efficacy of antibiotics, antiviral drugs, targeted cancer therapies, and other treatments. Accurately calculating protein–drug binding affinity changes upon mutations in the drug target is of high interest as this can yield a better understanding into how such mutations drive drug‐resistance, especially when the mutation in question does not directly interfere with binding of the drug. The main aim of this article is to provide an up‐to‐date reference on the computational tools that are available for the calculation of Gibbs energy (free energy) changes upon mutation, their strengths, and limitations. The methods that are discussed include free energy calculations (free energy perturbation, thermodynamic integration, multistate Bennett acceptance ratio), analysis of molecular dynamics simulations (linear interaction energy, molecular mechanics [MM]/Poisson–Boltzmann solvated area, and MM/generalized Born solvated area), and methods that involve quantum mechanical calculations (including QM/MM). The possibility to use machine learning is also introduced. Given that the benefit of accurately calculating binding affinity changes upon mutation depends on comparing calculated values with experimental measurements, a brief survey on experimental methods and observables is provided. Examples of computational studies that go beyond calculating the Gibbs energy changes are given. Factors that need to be addressed by the computational chemist and potential pitfalls are discussed at length.This article is categorized under: Structure and Mechanism > Computational Biochemistry and Biophysics Molecular and Statistical Mechanics > Free Energy Methods Molecular and Statistical Mechanics > Molecular Interactions

Coarse-Grained Modeling of Multiple Pathways in Conformational Transitions of Multi-Domain Proteins.

Large-scale conformational transitions in multi-domain proteins are often essential for their functions. To investigate the transitions, it is necessary to explore multiple potential pathways, which involve different intermediate structures. Here, we present a multi-basin (MB) coarse-grained (CG) structure-based Go̅ model for describing transitions in proteins with more than two moving domains. This model is an extension of our dual-basin Go̅ model in which system-dependent parameters are determined systematically using the multistate Bennett acceptance ratio method. In the MB Go̅ model for multi-domain proteins, we assume that intermediate structures may have partial inter-domain native contacts. This approach allows us to search multiple transition pathways that involve distinct intermediate structures using the CG molecular dynamics (MD) simulations. We apply this scheme to an enzyme, adenylate kinase (AdK), which has three major domains and can move along two different pathways. Using the optimized mixing parameters for each pathway, AdK shows frequent transitions between the Open, Closed, and the intermediate basins and samples a wide variety of conformations within each basin. The explored multiple transition pathways could be compared with experimental data and examined in more detail by atomistic MD simulations.

Extension of the Variational Free Energy Profile and Multistate Bennett Acceptance Ratio Methods for High-Dimensional Potential of Mean Force Profile Analysis.

We redevelop the variational free energy profile (vFEP) method using a cardinal B-spline basis to extend the method for analyzing free energy surfaces (FESs) involving three or more reaction coordinates. We also implemented software for evaluating high-dimensional profiles based on the multistate Bennett acceptance ratio (MBAR) method which constructs an unbiased probability density from global reweighting of the observed samples. The MBAR method takes advantage of a fast algorithm for solving the unbinned weighted histogram (UWHAM)/MBAR equations which replaces the solution of simultaneous equations with a nonlinear optimization of a convex function. We make use of cardinal B-splines and multiquadric radial basis functions to obtain smooth, differentiable MBAR profiles in arbitrary high dimensions. The cardinal B-spline vFEP and MBAR methods are compared using three example systems that examine 1D, 2D, and 3D profiles. Both methods are found to be useful and produce nearly indistinguishable results. The vFEP method is found to be 150 times faster than MBAR when applied to periodic 2D profiles, but the MBAR method is 4.5 times faster than vFEP when evaluating unbounded 3D profiles. In agreement with previous comparisons, we find the vFEP method produces superior FESs when the overlap between umbrella window simulations decreases. Finally, the associative reaction mechanism of hammerhead ribozyme is characterized using 3D, 4D, and 6D profiles, and the higher-dimensional profiles are found to have smaller reaction barriers by as much as 1.5 kcal/mol. The methods presented here have been implemented into the FE-ToolKit software package along with new methods for network-wide free energy analysis in drug discovery.

Variational Method for Networkwide Analysis of Relative Ligand Binding Free Energies with Loop Closure and Experimental Constraints.

We describe an efficient method for the simultaneous solution of all free energies within a relative binding free-energy (RBFE) network with cycle closure and experimental/reference constraint conditions using Bennett Acceptance Ratio (BAR) and Multistate BAR (MBAR) analysis. Rather than solving the BAR or MBAR equations for each transformation independently, the simultaneous solution of all transformations are obtained by performing a constrained minimization of a global objective function. The nonlinear optimization of the objective function is subjected to affine linear constraints that couple the free energies between the network edges. The constraints are used to enforce the closure of thermodynamic cycles within the RBFE network, and to enforce an additional set of linear constraint conditions demonstrated here to be subsets of (1 or 2) experimental values. We describe details of the practical implementation of the network BAR/MBAR procedure, including use of generalized coordinates in the minimization of the free-energy objective function, propagation of bootstrap errors from those coordinates, and performance and memory optimization. In some cases it is found that use of restraints in the optimization is more practical than use of generalized coordinates for enforcing constraint conditions. The fast BARnet and MBARnet methods are used to analyze the RBFEs of six prototypical protein-ligand systems, and it is shown that enforcement of cycle closure conditions reduces the error in the predictions only modestly, and further reduction in errors can be achieved when one or two experimental RBFEs are included in the optimization procedure. These methods have been implemented into FE-ToolKit, a new free-energy analysis toolkit. The BARnet/MBARnet framework presented here opens the door to new, more efficient and robust free-energy analysis with enhanced predictive capability for drug discovery applications.

Study of the effectiveness of various cannabinoid receptor 1 (CB1) agonists using molecular docking and molecular dynamics modeling

The binding of a series of small organic molecules, acting as agonists of the cannabinoid receptor CB1, was investigated by means of three methods of computational chemistry. Binding modes were predicted by means of molecular docking, and binding free energy was estimated via docking, molecular-mechanics Poisson-Boltzmann surface area method, and multistate Bennett acceptance ratio. No evident correlation was observed for the molecules between the experimental characteristics of affinity and three computed binding free energy estimates. The reasons for the discrepancy were discussed.

Building a Macro-mixing Dual-basin Go Model using the Multistate Bennett Acceptance Ratio

Molecular dynamics (MD) simulations of biomolecules are widely used to investigate conformational dynamics and structure-function relationship. All-atom (AA) models provide the most accurate description of the underlying dynamics. However at present, even with the fastest computers, the time-scale obtainable with an atomistic simulation is up to a millisecond, whereas biologically relevant motions occur at the time scale of milliseconds to seconds. To overcome this, coarse-grained (CG) modeling can be utilized. The use of CG models reduces the computational time by several orders of magnitude, allowing access to biologically relevant time-scales. In this work, we use the dual-basin Gō-model, which is a structural-based CG model, for simulating conformational transitions between two known structures of a protein. The dual-basin model is formed by mixing two single-basin potentials and includes system-dependent parameters which determine the barrier height and the energetic offsets between the two basins. The determination of parameters however, is usually not straightforward and can be time consuming. Here, we have developed an efficient scheme to determine the mixing parameters using the Multistate Bennett Acceptance Ratio (MBAR) method after short simulations with a set of parameters. In the scheme, MBAR allows us to predict observables at various unsimulated conditions, which are useful to improve the mixing parameters in the next round of iterative simulations. The number of iterations that are necessary for obtaining the converged mixing parameters is significantly reduced in the scheme. We applied the scheme to several proteins, for showing the effectiveness in parameter determination. After obtaining the converged parameters, the proteins show frequent conformational transitions between open and closed states, providing the theoretical basis to investigate structure-dynamics-function relationships of the proteins.

Overcoming Orthogonal Barriers in Alchemical Free Energy Calculations: On the Relative Merits of λ-Variations, λ-Extrapolations, and Biasing.

Alchemical free energy calculations using conventional molecular dynamics and thermodynamic integration rely on simulations performed at fixed values of the coupling parameter λ. When multiple conformers in equilibrium are separated by high barriers in the space orthogonal to λ, proper convergence may require extremely long simulations. Four main strategies can be employed to address this orthogonal-sampling problem: (a) λ-variations, where λ can change along the simulations to circumvent barriers; (b) λ-extrapolations, where statistical information is transferred between λ-points; (c) specific biasing, where orthogonal barriers are reduced using a biasing potential designed specifically for the system; and (d) generic biasing, where orthogonal barriers are reduced using a generic approach. Here, we investigate the relative merits of the first three strategies considering two benchmark systems. The KXK system involves a mutation of the central residue in a tripeptide to a glycine and the XTP system involves a hydrogen-to-bromine mutation in the base of a nucleotide. Three sampling methods are compared, the latter two involving λ-variations: molecular dynamics simulations with fixed λ-points, Hamiltonian replica exchange, and the recently introduced conveyor belt method. Two free energy estimators are applied, the second one involving λ-extrapolations: thermodynamic integration with Simpson quadrature and the multistate Bennett acceptance ratio. Finally, three different seeding schemes are considered for the generation of the initial configurations. For both benchmark systems, λ-extrapolations are found to provide little gain, whereas λ-variations can significantly enhance the convergence. They are sufficient on their own if the orthogonal barriers are low in at least one state (e.g., the glycine state in KXK). However, if the orthogonal barriers are high over the entire λ-range (e.g., the XTP system), λ-variations are only effective when applied together with a specific biasing for introducing such a low-barrier state.

Prediction of octanol-water partition coefficients for the SAMPL6-$$\log P$$ molecules using molecular dynamics simulations with OPLS-AA, AMBER and CHARMM force fields

All-atom molecular dynamics simulations with stratified alchemical free energy calculations were used to predict the octanol-water partition coefficient [Formula: see text] of eleven small molecules as part of the SAMPL6-[Formula: see text] blind prediction challenge using four different force field parametrizations: standard OPLS-AA with transferable charges, OPLS-AA with non-transferable CM1A charges, AMBER/GAFF, and CHARMM/CGenFF. Octanol parameters for OPLS-AA, GAFF and CHARMM were validated by comparing the density as a function of temperature, the chemical potential, and the hydration free energy to experimental values. The partition coefficients were calculated from the solvation free energy for the compounds in water and pure ("dry") octanol or "wet" octanol with 27 mol% water dissolved. Absolute solvation free energies were computed by thermodynamic integration (TI) and the multistate Bennett acceptance ratio with uncorrelated samples from data generated by an established protocol using 5-ns windowed alchemical free energy perturbation (FEP) calculations with the Gromacs molecular dynamics package. Equilibration of sets of FEP simulations was quantified by a new measure of convergence based on the analysis of forward and time-reversed trajectories. The accuracy of the [Formula: see text] predictions was assessed by descriptive statistical measures such as the root mean square error (RMSE) of the data set compared to the experimental values. Discarding the first 1 ns of each 5-ns window as an equilibration phase had a large effect on the GAFF data, where it improved the RMSE by up to 0.8 log units, while the effect for other data sets was smaller or marginally worsened the agreement. Overall, CGenFF gave the best prediction with RMSE 1.2 log units, although for only eight molecules because the current CGenFF workflow for Gromacs does not generate files for certain halogen-containing compounds. Over all eleven compounds, GAFF gave an RMSE of 1.5. The effect of using a mixed water/octanol solvent slightly decreased the accuracy for CGenFF and GAFF and slightly increased it for OPLS-AA. The GAFF and OPLS-AA results displayed a systematic error where molecules were too hydrophobic whereas CGenFF appeared to be more balanced, at least on this small data set.

Replica-Exchange Umbrella Sampling Combined with Gaussian Accelerated Molecular Dynamics for Free-Energy Calculation of Biomolecules.

We have developed an enhanced conformational sampling method combining replica-exchange umbrella sampling (REUS) with Gaussian accelerated molecular dynamics (GaMD). REUS enhances the sampling along predefined reaction coordinates, while GaMD accelerates the conformational dynamics by adding a boost potential to the system energy. The method, which we call GaREUS (Gaussian accelerated replica-exchange umbrella sampling), enhances the sampling more efficiently than REUS or GaMD, while the computational resource for GaREUS is the same as that required for REUS. The two-step reweighting procedure using the multistate Bennett acceptance ratio method and the cumulant expansion for the exponential average is applied to the simulation trajectories for obtaining the unbiased free-energy landscapes. We apply GaREUS to the calculations of free-energy landscapes for three different cases: conformational equilibria of N-glycan, folding of chignolin, and conformational change of adenyl kinase. We show that GaREUS speeds up the convergences of free-energy calculations using the same amount of computational resources as REUS. The free-energy landscapes reweighted from the trajectories of GaREUS agree with previously reported ones. GaREUS is applicable to free-energy calculations of various biomolecular dynamics and functions with reasonable computational costs.